//========================================================================================
// Sequence3Bit.hh
// ---------------------------------------------------------------------------------------
// Saves DNA sequences in 3 bits/base format.   
//
//----------------------------------------------------------------------------------------
// $Id: //CBT-depot/CBT++/include/Sequence3Bit.hh#4 $
//----------------------------------------------------------------------------------------
//
//========================================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8.........9
//


#include "CompBioTools.hh"

#ifndef _3BITSEQUENCE
#define _3BITSEQUENCE


#define PRECOMPUTE_TABLE_SIZE 5000
#define SEQ_BUFFER_SIZE3 500000


#define NORXASCIIMASK 0x08  // Mask to tell N/X directly from ASCII
#define NORX3BITMASK 0x04   // Mask to tell N/X from base in 3 bit encoding. 
#define ASCII3BITMASK 0x0E  // Mask to mask out bits 1,2,3 from ASCII

#define LOW3BITSMASK 0x00000007



//------------------------------------------------------------------------------
// * Sequence3Bit
//------------------------------------------------------------------------------
// A sequence class that stores bases in 3-bits. This allows n and x to be saved as well as 
// a,c,t,g while coming close to minimal space.  Since the goal is both space and speed 
// efficiency, specialized functions to read in Sequence3Bit sequences, or Sequence3BitVector 
// sets of sequences from a FASTA format file are included here, rather than as part of the 
// FASTAFile class.  
//
// Encoding is coordinated with the 2-bit encoding as follows:
//
// 2-bit/base encoding extracted from ASCII: 
// 
//      Hex         Dec            Binary   Encoding Bits
// a = 61 (hex)    97 (dec)      0110 0001   0110 0(00)1   a -> 0
// c = 63 (hex)    99 (dec)      0110 0011   0110 0(01)1   c -> 1
// g = 67 (hex)    103 (dec)     0110 0111   0110 0(11)1   g -> 3
// t = 74 (hex)    116 (dec)     0111 0100   0111 0(10)0   t -> 2
//
// A = 41 (hex)    65 (dec)      0100 0001   0100 0(00)1   A -> 0
// C = 43 (hex)    67 (dec)      0100 0011   0100 0(01)1   C -> 1
// G = 47 (hex)    71 (dec)      0100 0111   0100 0(11)1   G -> 3 
// T = 54 (hex)    84 (dec)      0101 0100   0101 0(10)0   T -> 2 
//
// So, bits 0,2 effectively code for a,c,t,g
//
// 3-bit encoding to include N/X.  
//
// When we use 3-bit bases, would like encoding for n/x to be a 
// superset of the acgtACGT encoding.   Bit 3 is set for n/x where it 
// is not set for acgtACGT, so bit 3 is diagnostic and bits 1,2,3 can 
// be used as the 3-bit encoding. 
//
// n = 6e (hex)    110 (dec)     0110 1110   0110 (111)0   n -> 7
// x = 78 (hex)    120 (dec)     0111 1000   0111 (100)0   x -> 4
//
// N = 4e (hex)    78 (dec)      0100 1110   0100 (111)0   N -> 7
// X = 58 (hex)    88 (dec)      0101 1000   0101 (100)0   X -> 4
// 
// 
//class Sequence3Bit : public vector<int>{
class Sequence3Bit {
 public:
  vector<int> mData;
  string mName;
  string mDescription;
  int    mBaseSize;
  int    mMerIdx;
  int    mBaseIdx;

  int  bases(){
    return(mBaseSize);
  }

 int size() const {
    return(mBaseSize);
  }

  void clear(){
    mData.clear();
    mName = "";
    mDescription = "";
    mBaseSize = 0;
    mMerIdx = 0;
    mBaseIdx = 0;
  }

  Sequence3Bit(){
    mBaseSize=0;mMerIdx = 0;mBaseIdx=0;	 
  }

  Sequence3Bit(Sequence &Seq){
    mBaseSize = 0;
    mBaseIdx = 0;
    mMerIdx = 0;
    mData.clear();
    ConcatenateSequence(Seq);
    mName = Seq.mName;
    mDescription = Seq.mDescription;
  }



  void ConcatenateSequence(Sequence &Seq);
  void ConcatenateCString(const char *cstring);
  void ConcatenateCStringTrimLowQual(const char *cstring);

  void Concatenate(Sequence3Bit &rhs);

  // In an effort to improve speed, precompute some of the 
  // info we need to convert a base index into the appropriate 
  // word and 3-bit offset. As it happens, though, this  was
  // a wash for performance. 
  //static int WordNum[PRECOMPUTE_TABLE_SIZE];
  //static int shiftAmt[PRECOMPUTE_TABLE_SIZE];
  //static unsigned int Mask3Bit[PRECOMPUTE_TABLE_SIZE];

  static char seqbuf3[SEQ_BUFFER_SIZE3];


  // deprecated. 
  int NumBases();
  
  void Init();

  Sequence3Bit &operator+(Sequence3Bit &rhs){
    (*this).Concatenate(rhs);
    return(*this);
  }

  // Make sure you can print out a sequence just like you can a string.
  friend ostream& operator<<(ostream& os,const Sequence3Bit &s){
    for(int i = 0;i<s.size();i++){
      os<<s[i];
    }
    return(os);
  } 

bool IsBaseNX(char base)
{
  switch (base)
    {
      case 'A': 
      case 'C': 
      case 'G': 
      case 'T': 
      case 'a': 
      case 'c': 
      case 'g': 
      case 't': 
      case 'x':
      case 'X':
      case 'N':
      case 'n':
	return(true);
	break;
      default:
	return(false);
    }
}

  Sequence3Bit& operator+=(char rhs){
    // KJD temporary hack. 
    
    if (IsBaseNX(rhs)){
      char buf[2];
      buf[0] = rhs;
      buf[1] = '\0';
      ConcatenateCString(buf);
    }
    return(*this);
  }
 
  Sequence3Bit& operator+=(Sequence3Bit &a){
    (*this).Concatenate(a);
    return(*this);
  }

  Sequence3Bit& operator=(const char*pSeqStr){
    (*this).ConcatenateCString(pSeqStr);
    return(*this);
  }


  // KJD Note:  formerly this was operator()
  // Want to explicitly not define operator() so it will
  // generate compile-time errors and we can fix other code. 

  //----------------------------------------
  // * operator[]
  //----------------------------------------
  //
  unsigned int operator[](int  n) const {
    // KJD debugging test remove when done. 
    // if (n > mBaseSize){
    //	Bug(n);
    //}else{		
    register unsigned int WordNum = n/10;		

    register unsigned int shiftAmt = (n - (WordNum*10))*3;
    register unsigned int mask = 0x07 << shiftAmt;
    return( (mData[WordNum] & mask) >> shiftAmt );
    // }
  } 


  //----------------------------------------
  // * Bit2Base
  //----------------------------------------
  // a = 61 (hex)    97 (dec)      0110 0001   0110 0(00)1   a -> 0
  // c = 63 (hex)    99 (dec)      0110 0011   0110 0(01)1   c -> 1
  // g = 67 (hex)    103 (dec)     0110 0111   0110 0(11)1   g -> 3
  // t = 74 (hex)    116 (dec)     0111 0100   0111 0(10)0   t -> 2
  // A = 41 (hex)    65 (dec)      0100 0001   0100 0(00)1   A -> 0
  // C = 43 (hex)    67 (dec)      0100 0011   0100 0(01)1   C -> 1
  // G = 47 (hex)    71 (dec)      0100 0111   0100 0(11)1   G -> 3 
  // T = 54 (hex)    84 (dec)      0101 0100   0101 0(10)0   T -> 2 
  // n = 6e (hex)    110 (dec)     0110 1110   0110 (111)0   n -> 7
  // x = 78 (hex)    120 (dec)     0111 1000   0111 (100)0   x -> 4
  // N = 4e (hex)    78 (dec)      0100 1110   0100 (111)0   N -> 7
  // X = 58 (hex)    88 (dec)      0101 1000   0101 (100)0   X -> 4  
  inline char Bit2Base(unsigned int x){
    switch(x){
      case 0: return('A');break;
      case 1: return('C');break;
      case 2: return('T');break;
      case 3: return('G');break;
      case 4: return('X');break;
      case 7: return('N');break;
      default: return('?');break;
    };
  }

  void Print(ostream &fout);
  void Set(int baseIdx,int merIdx, char c);
  bool GetNextReadTrimLowQual(istream &fin,char *buffer,int bufSize);
  bool GetNextRead(istream &fin,char *buffer,int bufSize);

  void GetSubSequence(Sequence3Bit &s,int start,int end);

  bool OnlyBases(int start,int chkSize);
  bool ContainsNX(int start,int chkSize);

  void ReverseComplement(Sequence3Bit &RevSeq);
  void Reverse(Sequence3Bit &RevSeq);
  void Reverse(){
    Reverse(*this);
  }

  Sequence3Bit(Sequence3Bit &S, int start, int end){
    S.GetSubSequence((*this),start,end);
  }

  void Bug(int n);

};


class Sequence3BitVector : public vector<Sequence3Bit>{
 public:

  void Read(istream &fin){
    Read(fin,readbuf,READ_BUFFER_SIZE);
  }

  void Read(istream &fin,char*buffer,int bufsize){	 
	 
    buffer[0] = '1';
    do{
      resize(size()+1);
      (*this)[size()-1].GetNextRead(fin,buffer,bufsize);
      //cout<<"Read "<<(*this)[size()-1].mName<<endl;
    }while(!fin.eof());
  }
};


#endif
